home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / qmult_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.5 KB  |  72 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Pre-multiply by random orthogonal matrix.
  4.  
  5. // Syntax:    B = qmult ( A )
  6.  
  7. // Description:
  8.  
  9. //    B is Q*A where Q is a random real orthogonal matrix from the
  10. //    Haar distribution, of dimension the number of rows in
  11. //    A. Special case: if A is a scalar then QMULT(A) is the same as
  12.  
  13. //        qmult(eye(a))
  14.  
  15. //       Called by RANDSVD.
  16.  
  17. //       Reference:
  18. //       G.W. Stewart, The efficient generation of random
  19. //       orthogonal matrices with an application to condition estimators,
  20. //       SIAM J. Numer. Anal., 17 (1980), 403-409.
  21.  
  22. //    This file is a translation of qmult.m from version 2.0 of
  23. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  24. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  25.  
  26. //-------------------------------------------------------------------//
  27.  
  28. qmult = function ( A )
  29. {
  30.   local (A)
  31.  
  32.   n = A.nr; m = A.nc;
  33.  
  34.   //  Handle scalar A.
  35.   if (max(n,m) == 1)
  36.   {
  37.     n = A;
  38.     A = eye(n,n);
  39.   }
  40.  
  41.   d = zeros(n,n);
  42.  
  43.   for (k in n-1:1:-1)
  44.   {
  45.     // Generate random Householder transformation.
  46.     rand("normal", 0, 1);
  47.     x = rand(n-k+1,1);
  48.     s = norm(x, "2");
  49.     sgn = sign(x[1]) + (x[1]==0);    // Modification for sign(1)=1.
  50.     s = sgn*s;
  51.     d[k] = -sgn;
  52.     x[1] = x[1] + s;
  53.     beta = s*x[1];
  54.  
  55.     // Apply the transformation to A.
  56.     y = x'*A[k:n;];
  57.     A[k:n;] = A[k:n;] - x*(y/beta);
  58.   }
  59.  
  60.   // Tidy up signs.
  61.   for (i in 1:n-1)
  62.   {
  63.     A[i;] = d[i]*A[i;];
  64.   }
  65.  
  66.   A[n;] = A[n;]*sign(rand());
  67.   B = A;
  68.  
  69.   rand("default");
  70.   return B;
  71. };
  72.